#include "BSpline.hpp"
#include "matplotlibcpp.h"
#include <cstdlib>
#include <cmath>
#include <vector>
#include <iomanip>

using namespace std;

double f(double x){
    return 1/(1+x*x);
}

double df(double x){
    return -2*x/pow(1+x*x,2);
}

void draw(vector<double> xData, vector<vector<double> > yData, string file_name){
    plt::figure_size(1200, 780); 
    plt::plot(xData, yData[0], {{"label", "f"}}); 
    plt::plot(xData, yData[1], {{"label", "Cubic_BSpline"}});
    plt::title("plot of C");
    plt::save(file_name,300);
    plt::legend();
}

int main() {
    std::cout << "----- running C and D -----" << std::endl;

    // C 的内容
    // 定义变量储存差值曲线
    vector<double> yValue;
    vector<double> xValue;
    double stepsize = 0.05;
    double current = -5;
    while(current <= 5){
        xValue.push_back(current);
        current += stepsize;
    }
    // 定义三次B样条，进行插值计算
    vector<double> xData;
    vector<double> yData;
    for (int i=1; i<=11; i++) {
        xData.push_back(-6+i);
        yData.push_back(f(-6+i));
    }
    double l_val = df(-5), r_val = df(5);
    BSpline c_bspline(3,xData,yData,BSpline::complete,l_val,r_val);
    c_bspline.setupSpline();
    vector<double> y_of_f, y_of_b;
    for (double x : xValue) {
        y_of_b.push_back(c_bspline.interpolate(x));
        y_of_f.push_back(f(x));
    }
    vector<vector<double> > y;
    y.push_back(y_of_f);
    y.push_back(y_of_b);
    draw(xValue, y, "../figure/C_cubic_plot.png");

    // 定义二次样条，进行插值计算
    xData.clear();
    yData.clear();
    yData.push_back(f(-5));
    for (int i=1; i<=10; i++) {
        //xData.push_back(i-5.5);
        yData.push_back(f(i-5.5));
    }
    yData.push_back(f(5));
    BSpline q_bspline(2);
    q_bspline.cardinalQuadratic(-5, 5, yData);
    stepsize = 0.05;
    current = -4.5;
    xValue.clear();
    while(current <= 4.5){
        xValue.push_back(current);
        current += stepsize;
    }
    y_of_b.clear();
    y_of_f.clear();
    for (double x : xValue) {
        y_of_f.push_back(f(x));
        y_of_b.push_back(q_bspline.interpolate(x));
    }
    y[0] = y_of_f;
    y[1] = y_of_b;
    draw(xValue, y, "../figure/C_quadratic_plot.png");

    // D的内容
    vector<double> c_error;
    vector<double> q_error;
    vector<double> points = {-3.5,-3,-0.5,0,0.5,3,3.5};
    for (double x : points) {
        c_error.push_back(abs(f(x) - c_bspline.interpolate(x)));
        q_error.push_back(abs(f(x) - q_bspline.interpolate(x)));
    }
    cout << "error as below:" << endl;
    // 设置输出格式
    std::cout << std::setw(10) << "Points" << " | " 
              << std::setw(15) << "C_Error" << " | " 
              << std::setw(15) << "Q_Error" << std::endl;

    // 输出points和对应的两种error
    for (size_t i = 0; i < points.size(); ++i) {
        std::cout << std::setw(10) << std::fixed << std::setprecision(1) << points[i] << " | " 
                  << std::setw(15) << std::fixed << std::setprecision(4) << c_error[i] << " | " 
                  << std::setw(15) << std::fixed << std::setprecision(4) << q_error[i] << std::endl;
    }

    std::cout << "----- complete C D -----" << std::endl;
}


